package oneD.fem.method;

import oneD.fem.structure.Structure;
import inf.jlinalg.IGESolver;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVectorReadOnly;
import inf.jlinalg.MatrixFormat;
import inf.jlinalg.QuadraticMatrixInfo;
import inf.jlinalg.SolveFailedException;
import inf.jlinalg.SymmetricGeneralGESolver;

public class ModalAnalyzer {
	private double[] omega;
	private IMatrix modalVectors;

	public ModalAnalyzer() {
	}

	public void analyze(Structure struct) throws SolveFailedException {
		int NEQ = struct.enumerateDOFs();
		// solver initialization
		IGESolver solver = new SymmetricGeneralGESolver();
		QuadraticMatrixInfo aInfo = solver.getAInfo();
		QuadraticMatrixInfo bInfo = solver.getBInfo();
		IMatrix kGlobal = solver.getA();
		IMatrix mGlobal = solver.getB();

		aInfo.setSize(NEQ);
		bInfo.setSize(NEQ);
		solver.initialize();

		struct.assembleStiffnessMatrix(kGlobal, true);
		struct.assembleMassMatrix(mGlobal);

		System.out.println("Stiffness matrix");
		System.out.println(MatrixFormat.format(kGlobal));
		System.out.println("Mass matrix");
		System.out.println(MatrixFormat.format(mGlobal));

		solver.solve();

		IVectorReadOnly lambda = solver.getEigenvalues();

		// read out eigen-frequency from eigenvalues
		omega = new double[NEQ];
		for (int i = 0; i < NEQ; i++)
			omega[i] = Math.sqrt(lambda.get(i));

		// extract modal vectors
		modalVectors = (IMatrix) solver.getEigenvectors();
		this.normalizeColumn(modalVectors);
	}

	public void printEigenfrequency() {
		System.out.println("idx\teigenfrequency");
		for (int i = 0; i < omega.length; i++)
			System.out.println(i + "\t" + omega[i]);
	}

	public void printModalVectors() {
		System.out.println("idx\tmodal vector");
		for (int i = 0; i < omega.length; i++) {
			System.out.print(i + "\t");
			for (int j = 0; j < omega.length; j++)
				System.out.print(modalVectors.get(j, i) + ",");
			System.out.println();
		}
	}

	private void normalizeColumn(IMatrix A) {
		for (int j = 0; j < A.getColumnCount(); j++) {
			double normColumn = 0.0;
			for (int i = 0; i < A.getRowCount(); i++)
				normColumn += Math.pow(A.get(i, j), 2);
			normColumn = Math.sqrt(normColumn);
			for (int i = 0; i < A.getRowCount(); i++)
				A.set(i, j, A.get(i, j) / normColumn);
		}
	}
}
